set more off

******************************************************************************************************************************
**********************************Figure 3, Table S.4, Table S.5***************************************************
******************************************************************************************************************************
*point estimate

local N = 54648
local J = 10
set seed 600
do programs/flexiblelogit_point.do `J'  `N' 
gen r=.
save point.dta, replace

********************************************************************************************************************************************* 
*********************************************beta confidence interval, Figure 3 (A), Table S.4************************************ 
*********************************************************************************************************************************************

use bootstrap,clear

append using point
forvalues i=1/3 {
gen impliedbeta`i'hatt=impliedbeta`i' if r==.
egen impliedbeta`i'hat=sum(impliedbeta`i'hatt)
drop impliedbeta`i'hatt
}
drop if r==. 
forvalues i=1/3{
sum impliedbeta`i'
gen var`i' =r(Var)
replace var`i'=1/var`i'
}
gen sumvar=var1+var2+var3

forvalues i=1/3{
gen weight`i'=var`i'/sumvar

}

append using point
forvalues i=1/3{
replace weight`i'=weight`i'[_n-1] if r==.

}
gen waveragebeta=weight1*impliedbeta1+weight2*impliedbeta2+weight3*impliedbeta3 

keep r waveragebeta 
gen wbetahatt=waveragebeta if r==. 
egen wbetahat=sum(wbetahatt)
drop wbetahatt 

*z0 
drop if r==.
gen low=(waveragebeta<=wbetahat)

egen sumlow=sum(low)

gen z0=invnormal(sumlow/250)

gen p1=normal(z0+z0-invnormal(0.975))
gen p2=normal(z0+z0+invnormal(0.975))

_pctile waveragebeta, nq(1000)
local p1p=max(1,round(p1*1000))
local p2p=min(999,round(p2*1000))
gen lower=r(r`p1p')
gen upper=r(r`p2p')

keep wbetahat lower upper 
keep if _n==1
rename (wbetahat lower upper)(beta_point beta_lb beta_ub)
label var beta_point "flexible logit point estimate of beta"
label var beta_lb "flexible logit lower bound of confidence interval"
label var beta_ub "flexible logit upper bound of confidence interval"


************************************************************************************************************
*************************Difference in beta absolute value, Figure 3(B), Table S.5**************************** 
************************************************************************************************************

use bootstrap,clear
append using point 


forvalues i=1/3 {
gen impliedbeta`i'hatt=impliedbeta`i' if r==.
egen impliedbeta`i'hat=sum(impliedbeta`i'hatt)
drop impliedbeta`i'hatt
}
drop if r==. 
forvalues i=1/3{
sum impliedbeta`i'
gen var`i' =r(Var)
replace var`i'=1/var`i'
}
gen sumvar=var1+var2+var3

forvalues i=1/3{
gen weight`i'=var`i'/sumvar

}
append using point 
forvalues i=1/3{
replace weight`i'=weight`i'[_n-1] if r==.

}
gen waveragebeta=weight1*impliedbeta1+weight2*impliedbeta2+weight3*impliedbeta3 

replace waveragebeta=abs(waveragebeta)-abs(naivez) 

keep r waveragebeta 
gen wbetahatt=waveragebeta if r==. 
egen wbetahat=sum(wbetahatt)
drop wbetahatt 

*z0 

drop if r==.
gen low=(waveragebeta<=wbetahat)
egen sumlow=sum(low)
gen z0=invnormal(sumlow/250)
gen p1=normal(z0+z0-invnormal(0.975))
gen p2=normal(z0+z0+invnormal(0.975))
_pctile waveragebeta, nq(1000)
local p1p=max(1,round(p1*1000))
local p2p=min(999,round(p2*1000))
gen lower=r(r`p1p')
gen upper=r(r`p2p')

keep wbetahat lower upper 
keep if _n==1
rename (wbetahat lower upper)(diffbeta_point diffbeta_lb diffbeta_ub)
label var diffbeta_point "point estimate of the difference between the absolute values of beta estimates from flexible logit and standard logit"
label var diffbeta_lb "lower bound of confidence interval"
label var diffbeta_ub "upper bound of confidence interval"



